class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##Módulo 2: Análisis espacial y Percepción Remota satelital óptica ###Introducción al análisis masivo: Paralelización ráster José A. Lastra<br> <a href="http://github.com/JoseLastra"> Github: JoseLastra</a><br> <a href="mailto:jose.lastra@pucv.cl"> jose.lastra@pucv.cl</a> | <a href="mailto:jose.lastramunoz@wur.nl"> jose.lastramunoz@wur.nl</a><br> .large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Octubre, 2024</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ## Contenidos .pull-left[ - Paralelización: * LibrerÃa parallel * LibrerÃa foreach y doParallel - LibrerÃa terra: * app() * sapp() * lapply() ] .pull-right[ <img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/> © Allison Horst ] --- ## Recordando algunos conceptos -- - Problemas de cálculo ([Jones, 2017](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html)): * **CPU**: proceso toma mucho tiempo de CPU. * **Memoria**: Demanda demasiada memoria. * **I/O**: toma mucho tiempo la lectura/escritura desde el disco. * **Network**: La red toma demasiado tiempo en transferir. -- - La paralelización viene a "solucionar" los problemas asociados a la CPU, dada la presencia de múltiples núcleos en los pc's modernos. --- ## Paralelización: Conceptos <center><img src="data:image/png;base64,#sock_fork.png" width="900px"/></center> <center>Fuente: Benito, 2021</center> --- ## foreach y doParallel: Ejemplo práctico con rásters -- - Crearemos una lista de archivos **MOD13Q1** correspondientes al Ãndice NDVI desde el año 2000 hasta el 2021. -- - Ejecutaremos dos tareas simples: (1) cortar y (2) enmascarar datos a un área de estudio. -- - Usaremos los siguientes archivos: (1) **san_javier.gpkg**, (2) **MOD13Q1_dates_ts492.csv** y (3) **MOD13Q1_bands (carpeta)** ``` r ## Listando archivos bandas <- list.files(path = "MOD13Q1_bands/", pattern = "*.tif", full.names = T) ## archivo de fechas tabla <- read_csv("MOD13Q1_dates_ts492.csv") ## shape de San Javier shp <- vect("san_javier.gpkg") ``` --- ## foreach y doParallel: Ejemplo práctico con rásters -- - Preparación de clúster de procesamiento. ``` r ## nombres de salida archivos nombres <- paste0( "San_Javier_MOD13Q1/", tabla$dates, "_", tabla$ID, "_San_Javier.tif" ) ## Cluster set up nCluster <- detectCores(logical = F) - 1 # Número de núcleos fÃsicos cl <- makeCluster(nCluster, type = "PSOCK") # Configuración del cluster registerDoParallel(cl) # registrar para doParallel ``` --- ## foreach y doParallel: Ejemplo práctico con raster -- - Configuración de foreach y detención del clúster ``` r foreach(i = 1:length(nombres)) %dopar% { # importar librerÃas a los nodos library(raster) library(magrittr) r <- bandas[i] %>% raster() # crear banda en R m <- r %>% crop(shp) %>% mask(shp) # #cortar y enmascarar # guardar resultado writeRaster(m, filename = nombres[i], overwrite = T, datatype = "INT2S" ) } stopCluster(cl) ``` -- - ¿Qué falló? --- ## foreach y doParallel: Ejemplo práctico con raster <center><img src="data:image/png;base64,#error.png" width="950px"/></center> -- - Esto lo podemos solucionar de varias maneras: - Utilizar __sf__ para leer nuestro vector: `shp <- read_sf('archivo.gpkg')` - Pasar la ruta del archivo y leerlo en cada procesamiento: `ruta <- 'san_javier.gpkg'` y `shp <- vect(ruta)` dentro del __foreach__. - Empaquetar nuestro SpatVector `shp <- vect('archivo.gpkg') %>% wrap()` y `vect(shp)` ó `unwrap(shp)` en nuestro __foreach__ --- <center><img src="data:image/png;base64,#cpu.png" width="600px"/></center> <center>LabGRS, 2024</center> --- ## terra y la paralelización -- - Los ejemplos de paralelización vistos anteriormente, donde se envÃa información a cada nodo, no son compatibles con los objetos de la librerÃa **terra** ([ver más Hijmans et al., 2024 pág. 8](https://cran.r-project.org/web/packages/terra/terra.pdf)) -- - terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo **app()** la más común y simple de emplear. -- - Podemos ajustar este comportamiento; pero algunas cosas no pueden ser trabajadas de esta forma. -- - Por esta razón, terra provee funciones que permiten realizar operaciones empleando múltiples núcleos. Siendo __app()__ la más común y directa. ``` r ndvi_sj <- list.files( path = "San_Javier_MOD13Q1/", pattern = "*.tif", full.names = T ) %>% rast() ``` --- ## app() -- - Esta función nos permite operar en un **SpatRaster** multi o monobanda para aplicar una función normalmente de resumen (sum, min, max, median, etc.). También podemos operar en **SpatRasterDataset**. -- - **Importante**: dado la eficiencia computacional de **terra** muchas funciones resumen pueden ser aplicadas directamente con un funcionamiento óptimo. -- - Sin embargo, para archivos muy grandes puede ser útil la aplicación de **app()** ``` r mean_noApp <- ndvi_sj %>% mean(na.rm = T) # simple mean_app <- app(ndvi_sj, fun = mean, na.rm = T, cores = 5) # 5 cores ``` -- - **app()** opera los archivos de la misma forma que **apply()**: donde cada banda es una columna en la matriz y los pÃxeles son las filas. --- <img src="data:image/png;base64,#DIPGEOPR_02_7_files/figure-html/unnamed-chunk-7-1.png" width="100%" /> --- ## sapp() -- - Si lo que queremos es aplicar una función a cada banda dentro de nuestro **SpatRaster**, podemos emplear *sapp()*. ``` r ndvi <- list.files( path = "MOD13Q1_bands/", pattern = "*.tif", full.names = T ) %>% rast() shp_v <- vect('san_javier.gpkg') # versión simple ndvi_sj <- ndvi %>% crop(shp_v) %>% mask(shp_v) ``` --- ## sapp() ``` r # versión con sapp() ## Creando la función crop_mask <- function(x, ...) { r <- x %>% crop(shp_v) %>% mask(shp_v) r } ## aplicación ndvi_sj_sapp <- sapp(ndvi_sj, fun = crop_mask) ``` --- <img src="data:image/png;base64,#DIPGEOPR_02_7_files/figure-html/unnamed-chunk-12-1.png" width="100%" /> --- ## Observaciones -- - **sapp()** invoca a **sapply()** por detrás y aún no se ha implementado la paralelización para esta alternativa (**Información actualizada al 14 de Octubre de 2024**). -- - Tal cuál lo indica la documentación: *"(...) La paralelización es soportada, pero a menudo no ayuda, e incluso puede ser más lenta."* -- - Justificaciones para paralelizar con terra: * Disponemos de más de 8 núcleos para procesar. * Tenemos un archivo muy grande que trabajar. * O tenemos una función muy compleja (lenta) para aplicar sobre nuestros archivos. * Muchas veces el uso de **lapply()** puede ser suficiente. --- ## lapply() con terra -- - Tomemos nuestra lista de archivos originales y la función creada anteriormente para cortar y enmascarar. ``` r ndvi_list <- list.files( path = "MOD13Q1_bands/", pattern = "*.tif", full.names = T ) # adaptando la función para lapply crop_mask <- function(x, ...) { x <- rast(x) # SpatRaster r <- x %>% crop(shp_v) %>% mask(shp_v) # aplicación de operación rm(x) # liberar memoria r # entregar objeto } ndvi_crop <- lapply(ndvi_list, FUN = crop_mask) ``` --- ## Recomendaciones -- - Siempre es mejor probar primero el código en un área pequeña antes de ejecutar el procesamiento de gran escala. -- - Podemos usar un área de prueba transformada en un tibble e ir porbando pixel a pixel la función, facilitando la búsqueda de errores. -- - También podemos emplear otras estructuras de código como transformar nuestros objetos en listas y data frames para paralelizar empleando __future_lapply()__ --- ## Funciones complejas en app() -- - El principal pro de __terra__ y `app()` tiene que ver con la potencia y eficiencia computacional. -- - Especialmente, si lo comparamos con su antecesor de librerÃa __raster__: `calc()` y `clusterR()` -- - Calculemos la tendencia de nuestro NDVI para el SpatRaster __ndvi_crop__ ``` r ## Creación de función para evaluar tendencia ## definiendo tiempo time <- 1:nlyr(ndvi_crop) calc_trend <- function(x, time) { if(all(is.na(x))){ # Si la serie es solo NA's devolvemos NA return(rep(NA,2)) } else { m <- lm(x ~ time)$coefficients[2] # extraccion pendiente p <- summary(lm(x ~ time))$coefficients[2,4] # extraccion pvalue out_vals <- c(m,p) # creacion de salida names(out_vals) <- c('trend', 'pval') # nombres para identificar variables return(out_vals) } } ``` --- --- ## Funciones complejas en app() ``` r ## aplicacion de funcion al archivo spatRaster NDVI_trend <- ndvi_crop %>% app(fun = calc_trend, time = time, cores = 6) ``` <center><img src="data:image/png;base64,#trend_pval.png" width="600px"/></center> --- ## BibliografÃa complementaria - Bosco, J. (2018). "R Para Principantes". [Ver](https://bookdown.org/jboscomendoza/r-principiantes4/) - Cao, R. y Fernández, R. (2022). "Introducción al procesamiento en paralelo en R". En: "Técnicas de Remuestreo". [Ver](https://rubenfcasal.github.io/book_remuestreo/intro-hpc.html) - Gillespie, C., & Lovelace, R. (2021). "Efficient R programming: a practical guide to smarter programming." O'Reilly Media, Inc. [Ver](https://csgillespie.github.io/efficientR/) - Hijmans, R. J., Bivand, R., Forner, K., Ooms, J., Pebesma, E., & Sumner, M. D. (2022). Package ‘terra’. [Ver](https://brieger.esalq.usp.br/CRAN/web/packages/terra/terra.pdf) - Jones, M. (2017). "Quick Intro to Parallel Computing in R". [Ver](https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html) - McCallum, E., & Weston, S. (2011). "Parallel R". O'Reilly Media, Inc. - Orellana, J. (2018). "HPC con R para investigadores". [Ver](https://bookdown.org/content/1498/) --- class: middle 